Predicting system trajectories toward critical transitions

ABSTRACT

Described is a system for predicting system trajectories toward critical transitions. The system transforms a set of multivariate time series of observables of a complex system into a set of symbolic multivariate time series. Then pair-wise time series of a transfer entropy (TE) measure are determined, wherein the TE measure quantifies the amount of information transfer from a source to a destination in the complex system. An associative transfer entropy (ATE) measure is determined which decomposes the pair-wise time series of TE to associative states of asymmetric, directional information flows, wherein the ATE measure is comprised of an ATE+ influence class and a ATE− influence class. The system estimates ATE+, TE, and ATE− trajectories over time, and at least one of the ATE+, TE, and ATE− trajectories is used to predict a critical transition in the complex system.

CROSS-REFERENCE TO RELATED APPLICATIONS

This is a Continuation-in-Part application of U.S. Non-Provisionalapplication Ser. No. 13/904,945, filed in the United States on May 29,2013, entitled, “Detection and Identification of Directional InfluencesUsing Dynamic Spectrum.” This is also a Non-Provisional patentapplication of U.S. Provisional Application No. 61/784,167, filed in theUnited States on Mar. 14, 2013, entitled, “Predict System TrajectoriesToward Critical Transitions.”

FIELD OF INVENTION

The present invention relates to system for automated detection ofcritical transitions and, more particularly, to a system for automateddetection of critical transitions using predicted system trajectories.

BACKGROUND OF INVENTION

One central characteristic of complex social-technological systems istheir self-organized emerging interactions which provide the benefits ofexchanging information and resources effectively, yet at the same timeincreases the risk and pace of spreading attacks or failures. A smallperturbation on a complex system operating in a high-risk unstableregion can induce a critical transition that leads to unstoppablecatastrophic failures.

There are a few recent works of early warning signals forheterogeneously networked dynamical systems. Kwon and Yang (see the Listof Incorporated Cited Literature References, Literature Reference No. 1)use transfer entropy to analyze financial market data. They calculatepair-wise transfer entropy to form a transfer entropy matrix anddemonstrate the asymmetrical influence from mature markets to emergingmarkets using transfer entropy. Their method estimates global transferentropy, not local transfer entropy changing in time; therefore, theirmethod cannot handle dynamics and structure changes in directionalinfluence.

Staniek and Lehnertz (see Literature Reference No. 4) introduced arobust and computationally efficient method to estimate transfer entropyusing a symbolization technique. Their demonstration of symbolictransfer entropy (STE) on brain electrical activity data shows STE isable to detect the asymmetric dependences and identify the hemispherecontaining the epileptic focus without observing actual seizureactivity. However, their emphasis is also not on identifying thedynamics of the complex system structure.

Scheffer et al. (see Literature Reference No. 5) surveysstate-of-the-art methods that originated in the ecological domain. Thesemethods focus only on homogeneous lattices with signals such asincreased temporal correlation, skewness, and spatial correlations ofthe population dynamics, thus are not able to deal with heterogeneouslynetworked dynamical systems.

Moon and Lu (see Literature References No. 2 and 3) developed a spectralearly warning signals (EWS) theory that detects the approaching ofcritical transitions and can estimate the system structure and networkconnectivity near critical transitions using the covariance matrix.Although spectral EWS quantifies how much entities of a system aremoving together, the symmetric nature of covariance spectrum does notpermit the analysis of directional influences among entities.

Thus, a continuing need exists for a system for automated detection ofcritical transitions in heterogeneously networked dynamical system usingan information dynamic spectrum framework that permits analysis ofdirectional influences among entities.

SUMMARY OF INVENTION

The present invention relates to system for predicting systemtrajectories toward critical transitions. The system comprises one ormore processors and a memory having instructions such that when theinstructions are executed, the one or more processors perform multipleoperations. A set of multivariate time series of observables of acomplex system is transformed into a set of symbolic multivariate timeseries. Pair-wise series of a transfer entropy (TE) measure aredetermined. An associative transfer entropy (ATE) measure is determined.The ATE measure is comprised of an ATE+ positive influence class and anATE− negative influence class. Then, ATE+, TE, and ATE− trajectories areestimated over time. Finally, a critical transition in the complexsystem is predicted using at least one of the ATE+, TE, and ATE−trajectories.

In another aspect, the system estimates the TE trajectory by the naturallogarithm (ln) with an unknown coefficient a and a constant c accordingto the following:

g(t)=a ln(t)+c,

where t represents time.

In another aspect, the system estimates the unknown coefficient a withvarious discrete time steps by taking the derivative of g(t)=a ln(t)+cto obtain

${{^{\prime}(t)} = \frac{a}{t}},$

then obtaining a discretized version of the unknown coefficient aaccording to the following:

$\frac{\Delta }{\Delta \; t} = {a{\frac{1}{t}.}}$

In another aspect, the system determines the unknown coefficient aaccording to the following:

-   -   for a fixed timestep Δt in a window [T_(start), T_(end)],        determining the following matrix equation:

${{\begin{bmatrix}\frac{1}{t_{1} + {\Delta \; t\text{/}2}} \\\frac{1}{t_{2} + {\Delta \; t\text{/}2}} \\\vdots \\\frac{1}{t_{k} + {\Delta \; t\text{/}2}}\end{bmatrix}a} = \begin{bmatrix}\frac{{\left( {t_{1} + {\Delta \; t}} \right)} - {\left( t_{1} \right)}}{\Delta \; t} \\\frac{{\left( {t_{2} + {\Delta \; t}} \right)} - {\left( t_{2} \right)}}{\Delta \; t} \\\vdots \\\frac{{\left( {t_{k} + {\Delta \; t}} \right)} - {\left( t_{k} \right)}}{\Delta \; t}\end{bmatrix}},$

where t₁, . . . , t_(k)+Δtε [T_(start), T_(end)]; and

-   -   determining an approximation of a_(Δt) for the fixed timestep Δt        according to the following:

${a_{\Delta \; t} = \left. {argmin}_{a}||{{\begin{bmatrix}\frac{1}{t_{1} + {\Delta \; t\text{/}2}} \\\frac{1}{t_{2} + {\Delta \; t\text{/}2}} \\\vdots \\\frac{1}{t_{k} + {\Delta \; t\text{/}2}}\end{bmatrix}a} - \begin{bmatrix}\frac{{\left( {t_{1} + {\Delta \; t}} \right)} - {\left( t_{1} \right)}}{\Delta \; t} \\\frac{{\left( {t_{2} + {\Delta \; t}} \right)} - {\left( t_{2} \right)}}{\Delta \; t} \\\vdots \\\frac{{\left( {t_{k} + {\Delta \; t}} \right)} - {\left( t_{k} \right)}}{\Delta \; t}\end{bmatrix}} \right.||_{2}},$

where argmin represents a minimization, t_(k) represents a k^(th) timepoint, and wherein the constant c corresponding to this timestep is

c _(Δt) =g(T _(end))−a _(Δt) ln(T _(end)).

In another aspect, the system determines a predicted value for a futuretime T_(end)+t_(d) according to the following:

-   -   F_(Δt)(t_(d))=a_(Δt) ln(T_(end)+t_(d))+c_(Δt), where t_(d) is        the desired time length for prediction.

BRIEF DESCRIPTION OF THE DRAWINGS

The objects, features and advantages of the present invention will beapparent from the following detailed descriptions of the various aspectsof the invention in conjunction with reference to the followingdrawings, where:

FIG. 1 is a block diagram depicting the components of a system forautomated detection of critical transitions system according to theprinciples of the present invention;

FIG. 2 is an illustration of a computer program product according to theprinciples of the present invention;

FIG. 3A illustrates a non-foster circuit in voltage over time accordingto prior art;

FIG. 3B illustrates transfer entropy (TE)/associative transfer entropy(ATE)+/ATE− curves over time according to prior an;

FIG. 4A illustrates curve fitting with a plot of the exponentialfunction according to the principles of the present invention;

FIG. 4B illustrates curve fitting with a plot of the logarithmicfunction according to the principles of the present invention;

FIG. 5 is a plot comparing an actual trajectory to a predictedtrajectory according to the principles of the present invention;

FIG. 6A illustrates a prediction cone produced according to theprinciples of the present invention and the predicted trajectory;

FIG. 6B illustrates a prediction cone produced according to theprinciples of the present invention and the predicted trajectory; and

FIG. 7 is a flow diagram depicting a system for automated detection ofcritical transitions according to the principles of the presentinvention.

DETAILED DESCRIPTION

The present invention relates to a system for automated detection ofcritical transitions and, more particularly, to a system for automateddetection of critical transitions using predicted system trajectories.

The following description is presented to enable one of ordinary skillin the art to make and use the invention and to incorporate it in thecontext of particular applications. Various modifications, as well as avariety of uses in different applications will be readily apparent tothose skilled in the art, and the general principles defined herein maybe applied to a wide range of aspects. Thus, the present invention isnot intended to be limited to the aspects presented, but is to beaccorded the widest scope consistent with the principles and novelfeatures disclosed herein.

In the following detailed description, numerous specific details are setforth in order to provide a more thorough understanding of the presentinvention. However, it will be apparent to one skilled in the art thatthe present invention may be practiced without necessarily being limitedto these specific details. In other instances, well-known structures anddevices are shown in block diagram form, rather than in detail, in orderto avoid obscuring the present invention.

The reader's attention is directed to all papers and documents which arefiled concurrently with this specification and which are open to publicinspection with this specification, and the contents of all such papersand documents are incorporated herein by reference. All the featuresdisclosed in this specification, (including any accompanying claims,abstract, and drawings) may be replaced by alternative features servingthe same, equivalent or similar purpose, unless expressly statedotherwise. Thus, unless expressly stated otherwise, each featuredisclosed is one example only of a generic series of equivalent orsimilar features.

Furthermore, any element in a claim that does not explicitly state“means for” performing a specified function, or “step for” performing aspecific function, is not to be interpreted as a “means” or “step”clause as specified in 35 U.S.C. Section 112, Paragraph 6. Inparticular, the use of “step of” or “act of” in the claims herein is notintended to invoke the provisions of 35 U.S.C. 112, Paragraph 6.

Before describing the invention in detail, first a list of citedreferences is provided. Next, a description of the various principalaspects of the present invention is provided. Subsequently, anintroduction provides the reader with a general understanding of thepresent invention. Finally, specific details of the present inventionare provided to give an understanding of the specific aspects.

(1) LIST OF INCORPORATED CITED LITERATURE REFERENCES

The following references are cited throughout this application. Forclarity and convenience, the references are listed herein as a centralresource for the reader. The following references are herebyincorporated by reference as though fully set forth herein. Thereferences are cited in the application by referring to thecorresponding literature reference number.

-   1. O. Kwon and J.-S. Yang, Information Flow between Stock Indices,    2008 EPL 82 68003.-   2. H. Moon and T.-C. Lu, Early Warning Signal of Complex Systems:    Network Spectrum and Critical Transitions, WIN 2010.-   3. H. Moon and T.-C. Lu, Network Catastrophe: Self-Organized    Patterns Reveal both the Instability and the Structure of Complex    Networks, preprint, 2012.-   4. M. Staniek and K. Lehnertz, Symbolic Transfer Entropy, Physical    Review Letters 100, 15801, 2008.-   5. M. Scheffer, J. Bascompte, W. A. Brock, V. Brovkin, S. R.    Carpenter, V. Dakos, H. Held, E. H. van Nes, M. Rietkerk, and G.    Sugihara, Early-warning signals for critical transitions. Nature,    461, 2009.-   6. T. Schreiber, Measuring Information Transfer, Phys. Rev, Lett.    85, 461, 2000.

(2) PRINCIPAL ASPECTS

The present invention has three “principal” aspects. The first is asystem for automated detection of critical transitions. The system istypically in the form of a computer system operating software or in theform of a “hard-coded” instruction set. This system may be incorporatedinto a wide variety of devices that provide different functionalities.The second principal aspect is a method, typically in the form ofsoftware, operated using a data processing system (computer). The thirdprincipal aspect is a computer program product. The computer programproduct generally represents computer-readable instructions stored on anon-transitory computer-readable medium such as an optical storagedevice, e.g., a compact disc (CD) or digital versatile disc (DVD), or amagnetic storage device such as a floppy disk or magnetic tape. Other,non-limiting examples of computer-readable media include hard disks,read-only memory (ROM), and flash-type memories. These aspects will bedescribed in more detail below.

A block diagram depicting an example of a system (i.e., computer system100) of the present invention is provided in FIG. 1. The computer system100 is configured to perform calculations, processes, operations, and/orfunctions associated with a program or algorithm. In one aspect, certainprocesses and steps discussed herein are realized as a series ofinstructions (e.g., software program) that reside within computerreadable memory units and are executed by one or more processors of thecomputer system 100. When executed, the instructions cause the computersystem 100 to perform specific actions and exhibit specific behavior,such as described herein.

The computer system 100 may include an address/data bus 102 that isconfigured to communicate information. Additionally, one or more dataprocessing units, such as a processor 104 (or processors), are coupledwith the address/data bus 102. The processor 104 is configured toprocess information and instructions. In an aspect, the processor 104 isa microprocessor. Alternatively, the processor 104 may be a differenttype of processor such as a parallel processor, or a field programmablegate array.

The computer system 100 is configured to utilize one or more datastorage units. The computer system 100 may include a volatile memoryunit 106 (e.g., random access memory (“RAM”), static RAM, dynamic RAM,etc.) coupled with the address/data bus 102, wherein a volatile memoryunit 106 is configured to store information and instructions for theprocessor 104. The computer system 100 further may include anon-volatile memory unit 108 (e.g., read-only memory (“ROM”),programmable ROM (“PROM”), erasable programmable ROM (“EPROM”),electrically erasable programmable ROM “EEPROM”), flash memory, etc.)coupled with the address/data bus 102, wherein the non-volatile memoryunit 108 is configured to store static information and instructions forthe processor 104. Alternatively, the computer system 100 may executeinstructions retrieved from an online data storage unit such as in“Cloud” computing. In an aspect, the computer system 100 also mayinclude one or more interfaces, such as an interface 110, coupled withthe address/data bus 102. The one or more interfaces are configured toenable the computer system 100 to interface with other electronicdevices and computer systems. The communication interfaces implementedby the one or more interfaces may include wireline (e.g., serial cables,modems, network adaptors, etc.) and/or wireless (e.g., wireless modems,wireless network adaptors, etc.) communication technology.

In one aspect, the computer system 100 may include an input device 112coupled with the address/data bus 102, wherein the input device 112 isconfigured to communicate information and command selections to theprocessor 100. In accordance with one aspect, the input device 112 is analphanumeric input device, such as a keyboard, that may includealphanumeric and/or function keys. Alternatively, the input device 112may be an input device other than an alphanumeric input device. In anaspect, the computer system 100 may include a cursor control device 114coupled with the address/data bus 102, wherein the cursor control device114 is configured to communicate user input information and/or commandselections to the processor 100. In an aspect, the cursor control device114 is implemented using a device such as a mouse, a track-ball, atrack-pad, an optical tracking device, or a touch screen. The foregoingnotwithstanding, in an aspect, the cursor control device 114 is directedand/or activated via input from the input device 112, such as inresponse to the use of special keys and key sequence commands associatedwith the input device 112. In an alternative aspect, the cursor controldevice 114 is configured to be directed or guided by voice commands.

In an aspect, the computer system 100 further may include one or moreoptional computer usable data storage devices, such as a storage device116, coupled with the address/data bus 102. The storage device 116 isconfigured to store information and/or computer executable instructions.In one aspect, the storage device 116 is a storage device such as amagnetic or optical disk drive (e.g., hard disk drive (“HDD”), floppydiskette, compact disk read only memory (“CD-ROM”), digital versatiledisk (“DVD”)). Pursuant to one aspect, a display device 118 is coupledwith the address/data bus 102, wherein the display device 118 isconfigured to display video and/or graphics. In an aspect, the displaydevice 118 may include a cathode ray tube (“CRT”), liquid crystaldisplay (“LCD”), field emission display (“FED”), plasma display, or anyother display device suitable for displaying video and/or graphic imagesand alphanumeric characters recognizable to a user.

The computer system 100 presented herein is an example computingenvironment in accordance with an aspect. However, the non-limitingexample of the computer system 100 is not strictly limited to being acomputer system. For example, an aspect provides that the computersystem 100 represents a type of data processing analysis that may beused in accordance with various aspects described herein. Moreover,other computing systems may also be implemented. Indeed, the spirit andscope of the present technology is not limited to any single dataprocessing environment. Thus, in an aspect, one or more operations ofvarious aspects of the present technology are controlled or implementedusing computer-executable instructions, such as program modules, beingexecuted by a computer. In one implementation, such program modulesinclude routines, programs, objects, components and/or data structuresthat are configured to perform particular tasks or implement particularabstract data types. In addition, an aspect provides that one or moreaspects of the present technology are implemented by utilizing one ormore distributed computing environments, such as where tasks areperformed by remote processing devices that are linked through acommunications network, or such as where various program modules arelocated in both local and remote computer-storage media includingmemory-storage devices.

An illustrative diagram of a computer program product (i.e., storagedevice) embodying the present invention is depicted in FIG. 2. Thecomputer program product is depicted as floppy disk 200 or an opticaldisk 202 such as a CD or DVD. However, as mentioned previously, thecomputer program product generally represents computer-readableinstructions stored on any compatible non-transitory computer-readablemedium. The term “instructions” as used with respect to this inventiongenerally indicates a set of operations to be performed on a computer,and may represent pieces of a whole program or individual, separable,software modules. Non-limiting examples of “instruction” includecomputer program code (source or object code) and “hard-coded”electronics (i.e. computer operations coded into a computer chip). The“instruction” is stored on any non-transitory computer-readable medium,such as in the memory of a computer or on a floppy disk, a CD-ROM, and aflash drive. In either event, the instructions are encoded on anon-transitory computer-readable medium.

(3) INTRODUCTION

One central characteristic of complex social-technological systems istheir self-organized emerging interactions which provide the benefits ofexchanging information and resources effectively, yet at the same timeincrease the risk and pace of spreading attacks or failures. A smallperturbation on a complex system operating in a high-risk unstableregion can induce a critical transition that leads to unstoppablecatastrophic failures. In the system according to the principles of thepresent invention, catastrophic failures can be avoided by detectingearly warnings of critical transitions and predicting the likelihood ofsystem trajectories and the lead time to the critical points.

The invention described herein is built upon the invention described inU.S. application Ser. No. 13/904,945, entitled, “Detection andIdentification of Directional Influences Using Dynamic Spectrum”(hereinafter referred to as the '945 application), which is aninformation dynamic spectrum framework that goes beyond unidirectionaldiffusion dynamics to investigate heterogeneously networked dynamicalsystems with transient directional influence. The '945 application ishereby incorporated by reference as though fully set forth herein.

The '945 application describes the development of a novel AssociativeTransfer Entropy (ATE) measure which decomposes the pair-wisedirectional influence of transfer entropy to associative states ofasymmetric, directional information flows. Multivariate time series ofcomplex systems were transformed into the spectrum of Transfer EntropyMatrix (TEM) and Associative Transfer Entropy Matrix (ATEM) to captureinformation dynamics of the system. The novel spectral radius measuresof TEM and ATEM were developed to detect early warning signals ofsource-driven instability, and induce directional influence structure toidentify the source and reveal dynamics of directional influences. Theframework described in the '945 application enables one not only toaccurately detect, predict, and mitigate abrupt behavior changes incomplex systems, but also to understand, design, and build moresustainable and resilient complex systems.

The idea of associative transfer entropy is to decompose transferentropy (TE) by constraining associated states of two processes. Thesimplest case of ATE is two influence classes: one positive ATE+ and theother negative ATE−. In this case, ATE is able to distinguish twosituations where the amount of information transfer is the same but haveopposite effects: the source drives the destination the same directionin ATE−, and the source drives the destination the opposite direction inATE−. For dynamic data, meaning the amount of information flow is notconstant, the TE/ATE of time series in a local time window iscalculated, so that TE/ATE becomes a function of time. Suchdecomposition further enables one to investigate ATEM in the form ofATEM+ and ATEM−.

Given ATE+/TE cross-over signature as the early warning signals ofcritical transition, the present invention predicts (1) the likelihoodof system trajectories and (2) the lead time to critical points. It wasfirst observed that ATE+ and TE are maximized, at the same time thatATE− flattened, prior to the system switching to alternative stableregimes. Therefore, model-based statistical forecasting was applied toestimate the likelihood of ATE+/TE/ATE− trajectories. The lead time ispredicted by estimating the critical point based on the projectedflattening of ATE− trajectories. The prediction method according to theprinciples of the present invention will continuously output the updateof probabilistic cones as the system progresses. Determining thelikelihood and lead time of critical points will enable furtherdevelopments of a mitigation strategy to avoid critical transitions.Each of these aspects will be described in further detail below.

(4) SPECIFIC DETAILS OF THE INVENTION

Described is a system to address the need of automated detection ofemerging transitions, and identify sources of instability from timeseries of complex systems. The techniques according to the principles ofthe present invention can be used to detect directional influences apartfrom intrinsic dynamics so that one can derive encapsulated interactionmodules to predict behavior trajectories. The system according to theprinciples of the present invention will enable the advancement inanalyzing dynamics of complex networks.

Advantages of the system and method described herein include, but arenot limited to, (1) providing an effective measure for quantifyingassociative, asymmetric directional influences, whereas other methodsonly can handle either symmetric or directional influences; (2)providing an effective formulation for capturing system-wise directionalinfluences, whereas other methods can only handle pair-wise directionalinfluence, or system-wise symmetrical influence; (3) providing aneffective measure for detecting instability in systems with directionalinfluence dynamics, whereas other methods assume un-directionaldiffusion dynamics; (4) providing an effective method for characterizinginstability trends of aggregated directional influences, whereas othermethods can only detect instability with simple thresholding; and (5)providing an effective method to reveal sources and structures inchanging dynamics of directional influences, whereas other methods canonly handle stationary dynamics.

(4.1) Background: Information Dynamic Spectrum

The asymmetric measure, transfer entropy,

$\begin{matrix}{{T_{x_{j}\rightarrow x_{i}} = {\Sigma_{{({x_{i,{t + 1}},x_{i,t},x_{j,t}})} \in D}{p\left( {x_{i,{t + 1}},x_{i,t},x_{j,t}} \right)}\log \frac{p\left( {\left. x_{i,{t + 1}} \middle| x_{i,t} \right.,x_{j,t}} \right)}{p\left( x_{i,{t + 1}} \middle| x_{i,t} \right)}}},} & (1)\end{matrix}$

quantifies the amount of information transfer from source x_(j) todestination x_(i). It has been used to demonstrate the asymmetricalinfluence from mature markets to emerging markets by analyzing stockmarket indices, as described in Literature Reference Nos. 1 and 6. Yet,current state-of-the-art only considers average pair-wise informationtransfers and statistical analysis, without including the morefundamental object of an information dynamic spectrum (i.e.,time-varying information transfer) to explain the emergence ofdirectional influences.

In the previously developed information dynamic spectrum theorydescribed in the '945 application, the dynamics of a complex networkwere formulated as: dX=F(X)dt+σdW, where F:

^(m)

^(m) is a differentiable function that explains the dynamics of thesystem of m elements X(t)=[x₁(t) x₂(t) . . . x_(m)(t)]^(T) and isdefined by F(X(t))=[f₁(X(t)) f₂(X(t)) . . . f_(m)(X(t))]^(T). Thiscaptures both the dynamics of individual entities and the networkedinteractions between the entities. For example, element x_(i)(t) can bethe i^(th) stock index of the t^(th) trading day. The last term, σdW,corresponds to white noise scaled by σ. To reveal the unknown F(X), thisproblem is further simplified by linear approximation:F(X)≈F(X₀)+J(X₀)(X−X₀), for X near X₀, where

$J = \left\lbrack \frac{\partial f_{i}}{\partial x_{j}} \right\rbrack_{i,{j = 1},{\ldots \; m}}$

is the Jacobian matrix. Therefore, the local temporal dynamics can beobtained pair-wise. The measure associative transfer entropy (ATE) atstate D_(k), for i≠j is defined by:

$\begin{matrix}{{T_{x_{j}\rightarrow x_{i}}^{D_{k}} = {\Sigma_{{({x_{i,{t + 1}},x_{i,t},x_{j,t}})} \in D}{p\left( {x_{i,{t + 1}},x_{i,t},x_{j,t}} \right)}\log \frac{p\left( {\left. x_{i,{t + 1}} \middle| x_{i,t} \right.,x_{j,t}} \right)}{p\left( x_{i,{t + 1}} \middle| x_{i,t} \right)}}},} & (2)\end{matrix}$

where D_(k) is a subset of D that represents a certain associated statebetween x_(i) and x_(j). The purpose of ATE is to capture informationtransfer between two variables for a particular state association. Forexample, TE may be decomposed into two influence classes: one positiveand the other negative.

If the internal dynamics within each node is negligible comparing tointeractions among nodes, it is conjectured that

${\frac{\partial f_{i}}{\partial x_{j}} \sim {\left( {T_{x_{j}\rightarrow x_{i}}^{D_{1}},T_{x_{j}\rightarrow x_{i}}^{D_{2}},\ldots,T_{x_{j}\rightarrow x_{i}}^{D_{s}}} \right)}},$

for some linear function g. Let T^(D) ^(k) be an associative transferentropy matrix (ATEM):

$\begin{matrix}{T^{D_{k}} = {\begin{bmatrix}T_{x_{1}\rightarrow x_{1}}^{D_{k}} & T_{x_{2}\rightarrow x_{1}}^{D_{k}} & \ldots & T_{x_{m}\rightarrow x_{1}}^{D_{k}} \\T_{x_{1}\rightarrow x_{2}}^{D_{k}} & T_{x_{2}\rightarrow x_{2}}^{D_{k}} & \ldots & T_{x_{m}\rightarrow x_{2}}^{D_{k}} \\\; & \vdots & \; & \; \\T_{x_{1}\rightarrow x_{m}}^{D_{k}} & T_{x_{2}\rightarrow x_{m}}^{D_{k}} & \ldots & T_{x_{m}\rightarrow x_{m}}^{D_{k}}\end{bmatrix}.}} & (3)\end{matrix}$

To appropriately deal with dynamic data, meaning the amount ofinformation flow from one to another is not static, TE/ATE curves weregenerated as functions of time, calculated over sliding windows.

Since the TEM matrix is non-symmetric, its eigenvalues arecomplex-valued. The spectral radius of the TEM matrix is used to measurethe total amount of information flow of the entire network.

To efficiently estimate transfer entropy of continuous data, Staniek andLehnertz (see Literature Reference No. 4) proposed a method using asymbolization technique that is robust and computationally fast, whichhad been introduced with the concept of permutation entropy. In the '945application, the symbolization technique of Staniek and Lehnertz wasadapted to calculate ATE. Specifically, the continuous-valued timesequence {x(t)}_(t=1) ^(N) was symbolized by first ordering the valuesof {x(t), x(t+1), . . . x(t+m)} with {1, 2, . . . , m} and denoting{circumflex over (x)}(t)=associated permutation of order m the symbol ofx(t) at t. Then, the ATE of {x(t)}_(t=1) ^(N) is estimated bycalculating ATE of {{circumflex over (x)}(t)}_(t=1) ^(N).

FIGS. 3A and 3B show an ATE analysis of a non-foster network from theinvention described in the '945 application. The circuit is initiallyoperated in the stable region, where there is no oscillation. Then asmall perturbation is added. It is unknown whether the circuit willbecome oscillatory (unstable) or stay stable. An ATE analysis wasperformed to detect if the circuits would become synchronized(unstable). The plot in FIG. 3A is the circuit in voltage over time,where Ant1 represents Antenna 1, Ant2 represents Antenna 2, i1represents port 1, and i2 represents port 2. The plot in FIG. 3B showsthe TE/ATE+/ATE− curves over time. The curves are obtained from theabsolute sum of spectrum of TEM/ATEM+/ATEM−, respectively. It was foundthat the spectral radius of TEM/ATEM+/ATEM− also show similar outcomes,but since TEM is the sum of ATEM+ and ATEM−, the ATE+ curve will nevercross the TE curve. Thus, the absolute sum of the spectrum is moreinformative. As shown in FIG. 3B, the ATE+ curve crosses over the TEcurve near time point 800, indicating the increasing in synchronization.In addition, the ATE+ curve reaches its peak right before fullsynchronization, while the ATE− curve flattens because there is nonegative association.

(4.2) Probabilistic Cones for Trajectories Prediction

In the system according to the principles of the present invention,model-based statistical forecasting is used to estimate the TE/ATE+/ATE−trajectories. An advantage of analyzing the TE/ATE curves is that TE/ATEremoves the spikes of the raw data, as shown in the circuit data examplein FIG. 3A. It has been observed that the TE curve in this case can bewell approximated by the natural logarithm with an unknown coefficient aand a constant c according to the following:

g(t)=a ln(t)+c.  (4)

FIG. 4A shows curve fitting with a plot of the exponential function, andFIG. 4B illustrates curve fitting with a plot of the logarithmicfunction. The bold solid curves 400 are the TE curve from FIG. 3B. Thebold dashed curves 402 are the fitted curves in the intervaltε[700,1100]. The unbolded dashed curves 404 are the projections fromthe fitted curves. FIG. 4B shows that the logarithm is appropriate tofit the TE curve.

To estimate the coefficient a from TE/ATE± curves, instead of fittingthem deterministically with the logarithmic function in equation (4),the rate of change was estimated with various discrete time steps. Thisgenerates a prediction cone, which will be described in detail below.From the derivative of equation (4),

$\begin{matrix}{{{^{\prime}(t)} = \frac{a}{t}},} & (5)\end{matrix}$

the following discretized version for the unknown a was obtained:

$\begin{matrix}{\frac{\Delta }{\Delta \; t} = {a{\frac{1}{t}.}}} & (6)\end{matrix}$

For a fixed timestep Δt in a window [T_(start), T_(end)], the leastsquares was used to solve the unknown a. First, the following matrixequation was written:

$\begin{matrix}{{{\begin{bmatrix}\frac{1}{t_{1} + {\Delta \; t\text{/}2}} \\\frac{1}{t_{2} + {\Delta \; t\text{/}2}} \\\vdots \\\frac{1}{t_{k} + {\Delta \; t\text{/}2}}\end{bmatrix}a} = \begin{bmatrix}\frac{{\left( {t_{1} + {\Delta \; t}} \right)} - {\left( t_{1} \right)}}{\Delta \; t} \\\frac{{\left( {t_{2} + {\Delta \; t}} \right)} - {\left( t_{2} \right)}}{\Delta \; t} \\\vdots \\\frac{{\left( {t_{k} + {\Delta \; t}} \right)} - {\left( t_{k} \right)}}{\Delta \; t}\end{bmatrix}},} & (7)\end{matrix}$

where t₁, . . . , t_(k)+Δtε[T_(start), T_(end)], where t_(k) is thek^(th) time point. The approximation of α_(Δt) for the fixed timestep Δtwas then obtained according to the following:

$\begin{matrix}{a_{\Delta \; t} = \left. {argmin}_{a}||{{\begin{bmatrix}\frac{1}{t_{1} + {\Delta \; t\text{/}2}} \\\frac{1}{t_{2} + {\Delta \; t\text{/}2}} \\\vdots \\\frac{1}{t_{k} + {\Delta \; t\text{/}2}}\end{bmatrix}a} - \begin{bmatrix}\frac{{\left( {t_{1} + {\Delta \; t}} \right)} - {\left( t_{1} \right)}}{\Delta \; t} \\\frac{{\left( {t_{2} + {\Delta \; t}} \right)} - {\left( t_{2} \right)}}{\Delta \; t} \\\vdots \\\frac{{\left( {t_{k} + {\Delta \; t}} \right)} - {\left( t_{k} \right)}}{\Delta \; t}\end{bmatrix}}||{}_{2}. \right.} & (8)\end{matrix}$

The constant corresponding to this timestep is thenc_(Δt)=g(T_(end))−a_(Δt) ln(T_(end)). Therefore, the predicted value forthe future time T_(end)+t_(d) is

F _(Δt)(t _(d))=a _(Δt) ln(T _(end) +t _(d))+c _(Δt).  (9)

t_(d) is the desired time length for prediction.

Now that a method to estimate the constants c and a has beenestablished, a probabilistic prediction cone (or light cone) can begenerated at each time point as one varies Δt to obtain multipleestimates of c and a. Therefore, a probabilistic prediction cone at agiven time consists of a collection of natural logarithm curves startingfrom that point.

To estimate the error of this prediction for the immediate nexttimestep, the following was calculated:

error_(Δt)(i)=g(t _(i) +Δt)−[a _(Δt) ln(t _(i) +Δt)+c(i)],  (10)

where

c(i)=g(t _(i−1))−a _(Δt) ln(t _(i−1)).  (11)

Let E be the collection of all error_(Δt) over all timesteps Δt and letσ=standard deviation of E. Let G be the collection of all predictedvalues G_(Δt)(t_(d)) over all timesteps Δt and μ=mean(G). Therefore, the95% confidence interval for the future time T_(end)+t_(d) is:

$\begin{matrix}{{{CI} = \left\lbrack {{\mu - {1.96\frac{\sigma}{\sqrt{N}}}},{\mu + {1.96\frac{\sigma}{\sqrt{N}}}}} \right\rbrack},} & (12)\end{matrix}$

where N is the size of E.

FIG. 5 shows the 95% confidence interval of a predicted trajectory inunbolded solid lines 500. The predicted trajectory for TE is representedby the unbolded dashed line 502. The bold solid line 504 depicts theactual trajectory of TE. The prediction value and confidence intervalare generated from the prediction cones, described above. Two points ofthe actual trajectory were outside the 95% confidence interval, rightbefore the non-foster circuits were fully synched around time t=1500.These points are represented by diamonds 504 and indicate thepossibility of a trend toward a critical transition.

FIGS. 6A and 6B illustrate two time snapshots of the prediction cones600 produced according to the method described above and the predictedtrajectory value (represented by bold dashed curves 602). The predictioncones are used to generate the prediction values with confidenceintervals in FIG. 5. The solid bold curves 604 represent the actual TEtrajectory.

FIG. 7 is a flow diagram depicting the system for predicting systemtrajectories toward critical transitions according to the principles ofthe present invention. In a first step 700, a set of multivariate timeseries of observables of a complex system is transformed into a set ofsymbolic multivariate time series. In a second step 702, the systemdetermines pair-wise time series of a TE measure. In a third step 704,an ATE measure is determined. In a fourth step 706, the system estimatesATE+, TE, and ATE− trajectories over time. In a fifth step 708, acritical transition in the system is predicted using at least one of theATE+, TE, and ATE− trajectories.

In summary, the present invention is a system and method for analyzingtime series of complex systems to detect trends toward criticaltransitions, especially for complex systems with directional influencedynamics. The described system and method can be deployed as embeddeddecision support modules to provide early indication of criticaltransitions in complex systems, including, but not limited to, vehiclehealth management, prognostics for complex electronics, crisismanagement in supply chains, social unrests, financial marketinstability, and epileptic seizures in brain networks. The technologydescribed herein will result in accurate detection of upcoming criticaltransitions of system behaviors, accurate characterization ofinstability trends toward critical transitions, and accurateidentification of instability sources and propagation structure changes.Applications of the present invention included, but are not limited to,vehicle health management, complex electronics, smart power grid, socialnetwork analysis, social instability, supply-chain crisis management,cyber security, and epileptic seizure detection.

1. A system for predicting system trajectories toward criticaltransitions, the system comprising: one or more processors installed asan embedded decision support module in a vehicle entity and anon-transitory computer-readable medium having executable instructionsencoded thereon such that when executed, the one or more processorsperform operations of: transforming a set of multivariate time series ofobservables of a complex system into a set of symbolic multivariate timeseries, wherein the complex system is a heterogeneously networkeddynamical system of vehicle entities; determining a transfer entropy(TE) measure between two time series, wherein the TE measure quantifiesan amount of information transfer from a source vehicle entity to adestination vehicle entity in the complex system; determining anassociative transfer entropy (ATE) measure by decomposing the TE measureto associative states of asymmetric, directional information flows, theassociative states being an ATE+ positive influence class and an ATE−negative influence class, wherein an influence from the source vehicleentity on the destination vehicle entity is positively correlated in theATE+ positive influence class such that if a value of the source vehicleentity is increasing, a value of the destination vehicle entity isincreasing, and an influence from the source vehicle entity on thedestination vehicle entity is negatively correlated in the ATE− negativeinfluence class such that if a value of the source vehicle entity isincreasing, a value of the destination vehicle entity is decreasing;estimating ATE+, TE, and ATE− trajectories over time; and predicting acritical transition in the complex system using at least one of theATE+, TE, and ATE− trajectories for analysis of directional informationinfluences among vehicle entities of the complex system to avoid acatastrophic failure of the complex system.
 2. The system as set forthin claim 1, wherein vehicle electrical systems the one or moreprocessors further perform an operation of estimating the TE trajectoryby the natural logarithm (ln) with an unknown coefficient a and aconstant c according to the following:g(t)=a ln(t)+c, where g(t) represents the estimated TE trajectory, andwhere t represents time.
 3. The system as set forth in claim 2, whereinthe one or more processors further perform an operation of estimatingthe unknown coefficient a with various discrete time steps by taking thederivative of g(t)=a ln(t)+c to obtain${{^{\prime}(t)} = \frac{a}{t}},$ then obtaining a discretized versionof the unknown coefficient a according to the following:$\frac{\Delta }{\Delta \; t} = {a{\frac{1}{t}.}}$
 4. The system asset forth in claim 3, wherein the one or more processors further performan operation of determining the unknown coefficient a according to thefollowing: for a fixed timestep Δt in a window [T_(start), T_(end)],determining the following matrix equation: ${{\begin{bmatrix}\frac{1}{t_{1} + {\Delta \; t\text{/}2}} \\\frac{1}{t_{2} + {\Delta \; t\text{/}2}} \\\vdots \\\frac{1}{t_{k} + {\Delta \; t\text{/}2}}\end{bmatrix}a} = \begin{bmatrix}\frac{{\left( {t_{1} + {\Delta \; t}} \right)} - {\left( t_{1} \right)}}{\Delta \; t} \\\frac{{\left( {t_{2} + {\Delta \; t}} \right)} - {\left( t_{2} \right)}}{\Delta \; t} \\\vdots \\\frac{{\left( {t_{k} + {\Delta \; t}} \right)} - {\left( t_{k} \right)}}{\Delta \; t}\end{bmatrix}},$ where t₁, . . . , t_(k)+Δtε[T_(start), T_(end)]; anddetermining an approximation of a_(Δt) for the fixed timestep Δtaccording to the following:${a_{\Delta \; t} = \left. \underset{a}{argmin}||{{\begin{bmatrix}\frac{1}{t_{1} + {\Delta \; t\text{/}2}} \\\frac{1}{t_{2} + {\Delta \; t\text{/}2}} \\\vdots \\\frac{1}{t_{k} + {\Delta \; t\text{/}2}}\end{bmatrix}a} - \begin{bmatrix}\frac{{\left( {t_{1} + {\Delta \; t}} \right)} - {\left( t_{1} \right)}}{\Delta \; t} \\\frac{{\left( {t_{2} + {\Delta \; t}} \right)} - {\left( t_{2} \right)}}{\Delta \; t} \\\vdots \\\frac{{\left( {t_{k} + {\Delta \; t}} \right)} - {\left( t_{k} \right)}}{\Delta \; t}\end{bmatrix}} \right.||_{2}},$ where argmin_(a) represents aminimization, t_(k) represents a k^(th) timestep, and wherein theconstant c corresponding to this timestep isc _(Δt) =g(T _(end))−a _(Δt) ln(T _(end)).
 5. The system as set forth inclaim 4, wherein the one or more processors further perform an operationof determining a predicted value for a future time T_(end)+t_(d)according to the following:F _(Δt)(t _(d))=a _(Δt) ln(T _(end) +t _(d))+c _(Δt), where F_(Δt) is apredicted value for a future time, and where t_(d) is a desired timelength for prediction.
 6. A computer-implemented method for predictingsystem trajectories toward critical transitions, comprising an act of:causing one or more processors installed as an embedded decision supportmodule in a vehicle entity to execute instructions stored on anon-transitory memory such that upon execution, the one or moreprocessors perform operations of: transforming a set of multivariatetime series of observables of a complex system into a set of symbolicmultivariate time series, wherein the complex system is aheterogeneously networked dynamical system of vehicle entities;determining a transfer entropy (TE) measure between two time series,wherein the TE measure quantifies an amount of information transfer froma source vehicle entity to a destination vehicle entity in the complexsystem; determining an associative transfer entropy (ATE) measure bydecomposing the TE measure to associative states of asymmetric,directional information flows, the associative states being an ATE+positive influence class and an ATE− negative influence class, whereinan influence from the source vehicle entity on the destination vehicleentity is positively correlated in the ATE+ positive influence classsuch that if a value of the source vehicle entity is increasing, a valueof the destination vehicle entity is increasing, and an influence fromthe source vehicle entity on the destination vehicle entity isnegatively correlated in the ATE− negative influence class such that ifa value of the source vehicle entity is increasing, a value of thedestination vehicle entity is decreasing; estimating ATE+, TE, and ATE−trajectories over time; and predicting a critical transition in thecomplex system using at least one of the ATE+, TE, and ATE− trajectoriesfor analysis of directional information influences among vehicleentities of the complex system to avoid a catastrophic failure of thecomplex system.
 7. The method as set forth in claim 6, wherein the oneor more processors further perform an operation of estimating the TEtrajectory by the natural logarithm (ln) with an unknown coefficient aand a constant c according to the following:g(t)=a ln(t)+c, where g(t) represents the estimated TE trajectory, andwhere t represents time.
 8. The method as set forth in claim 7, whereinthe data processor further performs an operation of estimating theunknown coefficient a with various discrete time steps by taking thederivative of g(t)=a ln(t)+c to obtain${{^{\prime}(t)} = \frac{a}{t}},$ then obtaining a discretized versionof the unknown coefficient a according to the following:$\frac{\Delta }{\Delta \; t} = {a{\frac{1}{t}.}}$
 9. The method asset forth in claim 8, wherein the data processor further performs anoperation of determining the unknown coefficient a according to thefollowing: for a fixed timestep Δt in a window [T_(start), T_(end)],determining the following matrix equation: ${{\begin{bmatrix}\frac{1}{t_{1} + {\Delta \; t\text{/}2}} \\\frac{1}{t_{2} + {\Delta \; t\text{/}2}} \\\vdots \\\frac{1}{t_{k} + {\Delta \; t\text{/}2}}\end{bmatrix}a} = \begin{bmatrix}\frac{{\left( {t_{1} + {\Delta \; t}} \right)} - {\left( t_{1} \right)}}{\Delta \; t} \\\frac{{\left( {t_{2} + {\Delta \; t}} \right)} - {\left( t_{2} \right)}}{\Delta \; t} \\\vdots \\\frac{{\left( {t_{k} + {\Delta \; t}} \right)} - {\left( t_{k} \right)}}{\Delta \; t}\end{bmatrix}},$ where t₁, . . . , t_(k)+Δt ε[T_(start), T_(end)]; anddetermining an approximation of a_(Δt) for the fixed timestep Δtaccording to the following:${a_{\Delta \; t} = \left. \underset{a}{argmin}||{{\begin{bmatrix}\frac{1}{t_{1} + {\Delta \; t\text{/}2}} \\\frac{1}{t_{2} + {\Delta \; t\text{/}2}} \\\vdots \\\frac{1}{t_{k} + {\Delta \; t\text{/}2}}\end{bmatrix}a} - \begin{bmatrix}\frac{{\left( {t_{1} + {\Delta \; t}} \right)} - {\left( t_{1} \right)}}{\Delta \; t} \\\frac{{\left( {t_{2} + {\Delta \; t}} \right)} - {\left( t_{2} \right)}}{\Delta \; t} \\\vdots \\\frac{{\left( {t_{k} + {\Delta \; t}} \right)} - {\left( t_{k} \right)}}{\Delta \; t}\end{bmatrix}} \right.||_{2}},$ where argmin_(a) represents aminimization, t_(k) represents a k^(th) timestep, and wherein theconstant c corresponding to this timestep isc _(Δt) =g(T _(end))−a _(Δt) ln(T _(end)).
 10. The method as set forthin claim 9, wherein the data processor further performs an operation ofdetermining a predicted value for a future time T_(end)+t_(d) accordingto the following:F _(Δt)(t _(d))=a _(Δt) ln(T _(end) +t _(d))+c _(Δt), where F_(Δt) is apredicted value for a future time, and where t_(d) is a desired timelength for prediction.
 11. A computer program product for predictingsystem trajectories toward critical transitions, the computer programproduct comprising computer-readable instructions stored on anon-transitory computer-readable medium that are executable by acomputer having one or more processors installed as an embedded supportmodule in a vehicle entity for causing the one or more processors toperform operations of: transforming a set of multivariate time series ofobservables of a complex system into a set of symbolic multivariate timeseries, wherein the complex system is a heterogeneously networkeddynamical system of vehicle entities; determining a transfer entropy(TE) measure between two time series, wherein the TE measure quantifiesan amount of information transfer from a source vehicle entity to adestination vehicle entity in the complex system; determining anassociative transfer entropy (ATE) measure by decomposing the TE measureto associative states of asymmetric, directional information flows, theassociative states being an ATE+ positive influence class and an ATE−negative influence class, wherein an influence from the source vehicleentity on the destination vehicle entity is positively correlated in theATE+ positive influence class such that if a value of the source vehicleentity is increasing, a value of the destination vehicle entity isincreasing, and an influence from the source vehicle entity on thedestination vehicle entity is negatively correlated in the ATE− negativeinfluence class such that if a value of the source vehicle entity isincreasing, a value of the destination vehicle entity is decreasing;estimating ATE+, TE, and ATE− trajectories over time; and predicting acritical transition in the complex system using at least one of theATE+, TE, and ATE− trajectories for analysis of directional informationinfluences among vehicle entities of the complex system to avoid acatastrophic failure of the complex system.
 12. The computer programproduct as set forth in claim 11, further comprising instructions forcausing the processor to perform an operation of estimating the TEtrajectory by the natural logarithm (ln) with an unknown coefficient aand a constant c according to the following:g(t)=a ln(t)+c, where g(t) represents the estimated TE trajectory, andwhere t represents time.
 13. The computer program product as set forthin claim 12, further comprising instructions for causing the processorto perform an operation of estimating the unknown coefficient a withvarious discrete time steps by taking the derivative of g(t)=a ln(t)+cto obtain ${{g^{\prime}(t)} = \frac{a}{t}},$ then obtaining adiscretized version of the unknown coefficient a according to thefollowing: $\frac{\Delta \; g}{\Delta \; t} = {a{\frac{1}{t}.}}$14. The computer program product as set forth in claim 13, furthercomprising instructions for causing the processor to perform anoperation of determining the unknown coefficient a according to thefollowing: for a fixed timestep Δt in a window [T_(start), T_(end)],determining the following matrix equation: ${{\begin{bmatrix}\frac{1}{t_{1} + {\Delta \; t\text{/}2}} \\\frac{1}{t_{2} + {\Delta \; t\text{/}2}} \\\vdots \\\frac{1}{t_{k} + {\Delta \; t\text{/}2}}\end{bmatrix}a} = \begin{bmatrix}\frac{{g\left( {t_{1} + {\Delta \; t}} \right)} - {g\left( t_{1} \right)}}{\Delta \; t} \\\frac{{g\left( {t_{2} + {\Delta \; t}} \right)} - {g\left( t_{2} \right)}}{\Delta \; t} \\\vdots \\\frac{{g\left( {t_{k} + {\Delta \; t}} \right)} - {g\left( t_{k} \right)}}{\Delta \; t}\end{bmatrix}},$ where t₁, . . . , t_(k)+Δtε[T_(start), T_(end)]; anddetermining an approximation of a_(Δt) for the fixed timestep Δtaccording to the following:${a_{\Delta \; t} = {\arg\limits_{a}\min \mspace{11mu} {{{\begin{bmatrix}\frac{1}{t_{1} + {\Delta \; t\text{/}2}} \\\frac{1}{t_{2} + {\Delta \; t\text{/}2}} \\\vdots \\\frac{1}{t_{k} + {\Delta \; t\text{/}2}}\end{bmatrix}a} - \begin{bmatrix}\frac{{g\left( {t_{1} + {\Delta \; t}} \right)} - {g\left( t_{1} \right)}}{\Delta \; t} \\\frac{{g\left( {t_{2} + {\Delta \; t}} \right)} - {g\left( t_{2} \right)}}{\Delta \; t} \\\vdots \\\frac{{g\left( {t_{k} + {\Delta \; t}} \right)} - {g\left( t_{k} \right)}}{\Delta \; t}\end{bmatrix}}}_{2}}},$ where argmin_(a) represents a minimization,t_(k) represents a k^(th) time point, and wherein the constant ccorresponding to this timestep isc _(Δt) =g(T _(end))−a _(Δt) ln(T _(end)).
 15. The computer programproduct as set forth in claim 14, further comprising instructions forcausing the processor to perform an operation of determining a predictedvalue for a future time T_(end)+t_(d) according to the following:F _(Δt)(t _(d))=a _(Δt) ln(T _(end) +t _(d))+c _(Δt), where F_(Δt) is apredicted value for a future time, and where t_(d) is a desired timelength for prediction.